Verify and Extend
2018-03-07
Verifiy and Extend
Verify and extend this, Jesse:
by my calculations, letting \(\hat{L}^{(R)}_d\) denote the resubstitution error rate of the quadratic Bayes plug-in classifier using \(d\) dimensions when \(f_0 = MVN([0,0],I)\) and \(f_1 = MVN([1,0],I)\), conditioning on \(n_0=n_1=n/2\) with \(n=20\) yields \[\mathbb{P}[ \hat{L}^{(R)}_2 > \hat{L}^{(R)}_1 ] \approx 0.18\] and \[\mathbb{P}[ \hat{L}^{(R)}_2 = \hat{L}^{(R)}_1 ] \approx 0.29\] that is, the probability that the second (useless) feature degrades performance, even when testing on the training set, is substantial.
Emprical results:
500 Monticarlo replicates with \(f_0 = MVN([0,0,0,0], I)\) and \(f_1 = MVN([1,0,0,0], I)\)
mont <- 500
d <- 4
n <- 20
n0 <- n1 <- n/2
truth <- c(rep(0,n0), rep(1,n1))
mu1 <- rep(0,d)
mu2 <- c(1, rep(0,d-1))
Ls <-
foreach(i = 1:mont, .combine = 'rbind') %dopar% {
set.seed(i)
s0 <- rmvnorm(n0, mean = mu1, sigma = diag(1,d))
s1 <- rmvnorm(n1, mean = mu2, sigma = diag(1,d))
samp_dat <- as.matrix(rbind(s0,s1))
Lrd <-
Reduce('cbind',
lapply(1:d, function(x){
sum(truth != predict(qda(truth ~ samp_dat[, 1:x]))$class)/length(truth)
}))
colnames(Lrd) <- paste0("Lr", 1:d)
as.data.frame(Lrd)
}| measurment | value |
|---|---|
| \(\mathbb{P}\left[\hat{L}^{(R)}_2 > \hat{L}^{(R)}_1\right]\) | 0.178 |
| \(\mathbb{P}\left[\hat{L}^{(R)}_2 = \hat{L}^{(R)}_1\right]\) | 0.296 |
| \(\mathbb{P}\left[\hat{L}^{(R)}_3 > \hat{L}^{(R)}_2\right]\) | 0.14 |
| \(\mathbb{P}\left[\hat{L}^{(R)}_3 = \hat{L}^{(R)}_2\right]\) | 0.252 |
| \(\mathbb{P}\left[\hat{L}^{(R)}_4 > \hat{L}^{(R)}_3\right]\) | 0.108 |
| \(\mathbb{P}\left[\hat{L}^{(R)}_4 = \hat{L}^{(R)}_3\right]\) | 0.282 |
| \(\mathbb{P}\left[\hat{L}^{(R)}_3 > \hat{L}^{(R)}_1\right]\) | 0.074 |
| \(\mathbb{P}\left[\hat{L}^{(R)}_3 = \hat{L}^{(R)}_1\right]\) | 0.174 |
| \(\mathbb{P}\left[\hat{L}^{(R)}_4 > \hat{L}^{(R)}_1\right]\) | 0.036 |
| \(\mathbb{P}\left[\hat{L}^{(R)}_4 = \hat{L}^{(R)}_1\right]\) | 0.076 |
Box for Carey; jitter and violin for Joshua
The plots below show the observed paired difference \(L^{(R)}_{(1:d)} - L^{(R)}_{1}\) for \(d = (2,3,4)\)
Showing the “going up” thing.
Here we chose a sub-sample where the \(L^{(R)}_{(1:d)} > L^{(R)}_{1}\) for \(d = (2,\dots,9)\) 
mu1 <- 1/2
ncDim <- 25
montc <- 500
nc <- 100
nc0 <- nc1 <- nc/2
truC <- c(rep(0,nc0), rep(1, nc1))
sc <- list()
for(mi in 1:montc){
tmp1 <- c(rnorm(nc0, mean = mu1, sd = 1), rnorm(nc1, mean =-mu1, sd = 1))
tmp2 <- sapply(1:(ncDim-1), function(x) rcauchy(nc, location = 0, scale = 1))
sc[[mi]] <- as.matrix(cbind(tmp1, tmp2))
}
qList <-
foreach(i = 1:100, .combine = 'rbind') %:%
foreach(j = 1:10, .combine = 'rbind') %dopar% {
dat <- sc[[i]][, 1:j]
qdaR <- qda(truC ~ dat)
rcl <- as.numeric(predict(qdaR)$class) - 1
qdaD <- qda(truC ~ dat, CV = TRUE)
(Lr <- sum(rcl != truC)/length(truC))
(Ld <- sum((as.numeric(qdaD$class) - 1) != truC)/length(truC))
rbind(
data.frame(Type = "Lr", value = Lr, Run = i, Dim = j),
data.frame(Type = "Ld", value = Ld, Run = i, Dim = j))
}
#m1 <- melt(qList[, 1:2])
#m1 <- melt(qList)
p <- ggplot(qList, aes(x = Type, y = value)) + geom_boxplot(notch = TRUE, alpha = 0.75) +
geom_jitter(aes(color = Type), alpha = 0.5, size = 0.5) +
facet_grid(. ~ Dim) + ggtitle("1st normal, then Cauchy")
show(p)## Warning: Removed 500 rows containing non-finite values (stat_boxplot).## notch went outside hinges. Try setting notch=FALSE.
## notch went outside hinges. Try setting notch=FALSE.## Warning: Removed 500 rows containing missing values (geom_point).
pairs(sc[[1]][, 1:5], col = alpha(truC + 1,0.5), pch = 20, cex = 0.4, xlim = c(-2,2), ylim = c(-2,2))
mu1 <- 1/2
ncDim <- 25
montc <- 500
nc <- 100
nc0 <- nc1 <- nc/2
truC <- c(rep(0,nc0), rep(1, nc1))
sc <- list()
for(mi in 1:montc){
tmp1 <- c(rnorm(nc0, mean = mu1, sd = 1), rnorm(nc1, mean =-mu1, sd = 1))
tmp2 <- sapply(1:(ncDim-1), function(x) rnorm(nc, mean = 0, sd = sqrt(1/128)))
sc[[mi]] <- as.matrix(cbind(tmp1, tmp2))
}
aList <-
foreach(i = 1:montc, .combine = 'rbind') %:%
foreach(j = 1:ncDim, .combine = 'rbind') %dopar% {
dat <- sc[[i]][, 1:j]
qdaR <- qda(truC ~ dat)
rcl <- as.numeric(predict(qdaR)$class) - 1
qdaD <- qda(truC ~ dat, CV = TRUE)
(Lr <- sum(rcl != truC)/length(truC))
(Ld <- sum((as.numeric(qdaD$class) - 1) != truC)/length(truC))
rbind(
data.frame(Type = "Lr", value = Lr, Run = i, Dim = j),
data.frame(Type = "Ld", value = Ld, Run = i, Dim = j))
}
#m1 <- melt(qList[, 1:2])
#m1 <- melt(qList)
p2 <- ggplot(aList, aes(x = Type, y = value)) +
geom_jitter(aes(color = Type), alpha = 0.5, size = 0.25) +
geom_boxplot(notch = TRUE, alpha = 0.5) +
facet_grid(. ~ Dim) + ggtitle("1st normal (0.5,-0.5), rest normal(0,0)")
show(p2)## notch went outside hinges. Try setting notch=FALSE.
## notch went outside hinges. Try setting notch=FALSE.
pairs(sc[[1]][, 1:10], col = alpha(truC + 1,0.2), pch = 20, cex = 0.4, xlim = c(-2,2), ylim = c(-2,2))
plot(sc[[1]][,1])
mu1 <- 1/2
ncDim <- 500
montc <- 100
nc <- 5000
nc0 <- nc1 <- nc/2
truC <- c(rep(0,nc0), rep(1, nc1))
sc <- list()
set.seed(317)
for(mi in 1:montc){
tmp1 <- c(rnorm(nc0, mean = mu1, sd = 1), rnorm(nc1, mean =-mu1, sd = 1))
tmp2 <- sapply(1:(ncDim-1), function(x) rnorm(nc, mean = 0, sd = sqrt(1/128)))
sc[[mi]] <- as.matrix(cbind(tmp1, tmp2))
}
bList <-
foreach(i = 1:50, .combine = 'rbind') %:%
foreach(j = c(1,seq(50,ncDim, by = 50)), .combine = 'rbind') %dopar% {
dat <- sc[[i]][, 1:j]
qdaR <- qda(truC ~ dat)
rcl <- as.numeric(predict(qdaR)$class) - 1
qdaD <- qda(truC ~ dat, CV = TRUE)
(Lr <- sum(rcl != truC)/length(truC))
(Ld <- sum((as.numeric(qdaD$class) - 1) != truC)/length(truC))
rbind(
data.frame(Type = "Lr", value = Lr, Run = i, Dim = j),
data.frame(Type = "Ld", value = Ld, Run = i, Dim = j))
}
#m1 <- melt(qList[, 1:2])
#m1 <- melt(qList)
p3 <- ggplot(bList, aes(x = Type, y = value)) + geom_boxplot(notch = TRUE) +
geom_jitter(aes(color = Type), alpha = 0.5, size = 0.25) +
facet_grid(. ~ Dim) + ggtitle("1st normal (0.5,-0.5), rest normal(0,0)")show(p3)## notch went outside hinges. Try setting notch=FALSE.
pairs(sc[[1]][, 1:10], col = alpha(truC + 1,0.2), pch = 20, cex = 0.4, xlim = c(-2,2), ylim = c(-2,2))